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Abstract — For reconstruction of low-rank matrices from under- 
sampled measurements, we develop an iterative algorithm based 
on least-squares estimation. While the algorithm can be used 
for any low-rank matrix, it is also capable of exploiting a-priori 
knowledge of matrix structure. In particular, we consider linearly 
structured matrices, such as Hankel and Toeplitz, as well as 
positive semidefinite matrices. The performance of the algorithm, 
referred to as alternating least-squares (ALS), is evaluated by 
simulations and compared to the Cramer-Rao bounds. 

Index Terms — Low-rank matrix reconstruction, Cramer-Rao 
bound, least squares, structured matrices 



I. Introduction 

Low-rank matrices appear in various areas of signal pro- 
cessing and system identification (TJ. In recent times the 
problem of reconstructing such matrices from a vector of 
linear measurements embedded in noise has attracted a lot 
of attention Q. In this scenario the dimension of the mea- 
surement vector is lower than the number of matrix elements, 
and hence the problem consists of an underdetermined set of 
linear equations. This problem has a myriad of applications, 
for example, spectral imaging 0, wireless sensor networks 
(4), video error concealment 0. 

In the gamut of designing low-rank matrix reconstruction 
algorithms, most of the existing research work deals with a 
specific setup known as 'matrix completion' where the mea- 
surement vector consists of a subset of elements of the under- 
lying matrix. The algorithms can be separated into three broad 
categories: Convex-relaxation based El, Q, minimization on 
Grassmannian manifold of subspaces JS), 0, least-squares 
matrix fitting iflOl . ifTTl . While we note that a substantial effort 
is paid to the matrix completion problem, far less effort is 
devoted to a general setup consisting of an underdetermined 
system of linear equations. There are few attempts, such as 
£3, d, E) and HI. 

For the general underdetermined setup, we first consider the 
aspect of designing an algorithm that can deal with any low- 
rank matrix. In addition, we consider the aspect of exploiting 
a-priori knowledge of underlying matrix structure for further 
performance improvement. The motivation for using structured 
low-rank matrices is due to their frequent occurrence in 
various signal processing and system identification problems. 
We develop a new algorithm that addresses both aspects. 
Extending the approach of IfTTl . the new algorithm is devel- 
oped in an iterative framework based on least-squares (LS) 
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estimation. For investigating matrix structure, we consider 
linearly structured matrices, such as Hankel and Toeplitz, as 
well as positive semidefinite matrices. The performance of the 
algorithm is evaluated by simulations and then compared to 
the Cramer-Rao bounds (CRBs). The CRBs are presented for 
measurements in Gaussian noise. Specifically we derived the 
CRBs for Hankel and positive semidefinite matrices. 

Notation: The invertible vectorization and matrix construc- 
tion mappings are denoted vec(-) : C" xp — > C" pxl and 
mat n , p (-) : C npxl -> C nxp , respectively. X r = {X e 
C nx P . rank ( X ) = r } and X + = {X e C nxn : X t 0} 
denote the sets of rank r matrices and positive semidefinite 
(p.s.d.) matrices, respectively. Similarly, Xg = {X 6 c nx P : 
vec(X) = S0,0 € C 9 } denotes a subspace of linearly 
structured matrices specified by S 6 i^npxq^ wn j c h includes 
Hankel and Toeplitz matrices. A<g)B is the Kronecker product 
and (A, B) = tr (B* A) is the inner product. A T , A* and At 
denote the transpose, Hermitian transpose and Moore-Penrose 
pseudoinverse of A e C nxp , respectively while ||A||f is its 
Frobenius norm. ||x||w = Vx*Wx is the weighted norm. 

II. System model 

A. General underdetermined setup 

We consider a matrix X e X r , where r <C min(n,p) is 
assumed to be known. It is observed by an undersampled linear 
measurement process, 

y = ^(X) + neC mxl , (1) 

where m < np and the linear sensing operator A : C nxp — > 
C mxl can be written equivalently in forms, 



(X,Ai) 



(X,A ?1 



Avec(X). 



(2) 



The matrix A is assumed to be known and the measurement 
noise n is assumed to be zero-mean, with E[nn*] = C € 
C mxm given. Note that in matrix completion, {A^} is nothing 
but the set of element-selecting operators. 

B. Applications in signal processing and system identification 

Applications of the general underdetermined setup are il- 
lustrated in the following examples: 

1) Recovery of data matrices X € X r , compressed by some 
randomly chosen linear operation A into y. 

2) Reconstruction of covariance matrices Cov(z) = R S 
X r n X+ from a subset of second-order moments = 



2 



E[zjZ?], estimated with zero-mean error £y, 
Then (Q~|) can be applied, y 



(i,j)efi. 

«4(R) + e. In certain 
applications R G X r (1 n A5, e.g. Toeplitz or 
Persy mmetric structure. 

3) Recovery of distance matrices D G AV, where dji — 
dij = ||xj — XjH^f and W >- 0. A subset of distance 
measurements are observed in noise. 

4) Identification of a system matrix D e A" r . The system 
z t = Du ( G C mxl is sampled by a varying linear 
operator y f = A t z t + n t , where A t G C sxp . A special 



case is A t 



-Ht)< 



where the index £(t) varies the 



sampling of zt at each t. Given a set of input and output 
data, {yt,Ut}J =1 , where sT < np, the observations can 
be stacked by vectorization, 



u. 



vec(D) 



n. 



In certain scenarios D may also have linear structure. 

III. Alternating least-squares estimator 
A. General low-rank matrix reconstruction 

We begin by considering the case of reconstructing X G X r 
Using a weighted least-squares criterion, the estimator is 



X = argmin||y-^(X)||^. 



(3) 



When the measurement noise n is Gaussian, this estimator 
coincides with the maximum likelihood estimator. For brevity 
we assume spatially uncorrelated noise, C = u 2 I m , without 
loss of generality. Then minimizing the ^2-norm is equivalent 
to (01. For general C the observation model is pre-whitened 
by forming y = C _1 / 2 y and A = C _1 / 2 A. 

Since X G X r , we express X = LR where L G C nxr and 
R G C rxp . Then the square of the measurement residual can 
be written as 



J(L,R)4|| y -^(LR)|| 
= ||y- A(I p (§ 

= lly-A(R T 



L)vec(R)|| 2 
»I„)vec(L)| 



(4) 



The cost J(L, R) is minimized in an alternating fashion by 
the following steps: 

• minimizing R with a fixed L, 

• minimizing L with a fixed R. 

In the new algorithm, the alternating minimization is per- 
formed through iterations. Starting with an initial L, the 
iterations continue as long as the decreasing trend of J(L, R) 
is observed. Given L, the minimizer of R is computed by 
vec(R) = [A(rp®L)]Ty and similarly, given R, the minimizer 
of L is computed by vec(L) = [A(R T <E> I„)]Ty. Since 
the Kronecker products are sparse, A(I P ® L) G C mxrp 
and A(R T ® I n ) G C mxnr can be computed efficiently. 
Henceforth, we refer to the algorithm as alternating LS (ALS). 



B. Structured low-rank matrix reconstruction 

Next, we consider a structured low -rank matrix X G X r , and 
develop an ALS for a known matrix structure in Algorithm [T] 
In the algorithm, for each iteration, we approach the LS 
problem by first relaxing the structural constraint, and compute 
R with a fixed L. Then, to impose the structural constraint 
on R, the low-rank matrix estimate is projected onto the set 
of structured matrices by X = 'P(LR), similar to 'lift and 
project' lfl5l . R is subsequently modified as the least-squares 
solution of X, 



min 1 1 LR - 

R 



x|||. 



L is updated in the same fashion. Here we mention that the 
algorithm is initialized by L := U0S0 where [Uo,Xo,Vo] 
= svdtrunc(mat„,p(A*y), r) and svdtrunc(Z, r) denotes 
the singular value decomposition of Z G C nxp truncated to 
the rth singular value. 

Algorithm 1 : ALS with known matrix structure 
1 
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Input: y, A and r 

[U ,So,V ] = svdtrunc(mat„ iP (A*y),r) 
L:=U £ 

while J(L, R) decreases do 



:= mat.^ ([A(I P (g 
:=P(LR) 
:= L-tX 

= mat„.. r ([A(R T 
:=7>(LR) 
= XRt 
end while 
Output: X = LR 



R 

X 
R 
L : 

X 

L : 



L)] f y) 



i n )] f y) 



For linearly structured matrices the projection, X = 
"P;t s (LR), is computed by obtaining 9 = S^vec(LR) and 
setting X = mat„. p (S0). For p.s.d. matrices, X = "P;f + (LR) 
is computed by first performing an eigenvalue decomposition 
of a symmetrized matrix, which projects the matrix onto the 
set of Hermitian matrices, ±(LR + (LR)*) = VAV*. The 
decomposition is performed to the rth largest eigenvalue, then 
setting X = V r A r V*, where positive eigenvalues have been 
retained, projects the matrix onto X+. 

IV. Cramer-Rao bounds 
For relevant comparison, we use CRBs. In this section, we 
describe the CRB expressions. 

A. Unstructured matrix 

For real-valued X and Gaussian measurement noise n, y 
is distributed as A/"(Avec(X), C). The Cramer-Rao bound for 
unbiased low-rank matrix estimators without structure was de- 
rived in HQ (see also Q): E J/ | X [||X-X(y)|||] > CRB(X), 
where 

CRB(X) = tr ((P T A T C" 1 AP)- 1 ) . (5) 

The CRB holds when rank(AP) = r(n + p) — r 2 where 
P = [Vi ® U U ® V V ® Ui] . The submatrices are 
obtained from the left-singular vectors, U = [u 1; ...,u r ] 
and Ui = [u r+1 , . . . , u„], and right-singular vectors, V = 
[vi, v r ] and Vi = [v r+ i, . . . , v„], of X G X r . 
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B. Structured matrices 

1) Hankel and Toeplitz matrices: For certain linearly 
structured matrices, such as Hankel or Toeplitz, low-rank 
parametrizations SO = Sg(a) exist that enable the derivation 
of CRBs. E.g. when X$ is Hankel, X can be parameterized in 
the canonical controllable form 1161 as [X]^ = b T $ 4+ -' _2 e 1 , 
where b e K r , 



and ei is the first standard basis vector in K r . Here the 
parameters are a = [a T b T ] , where a = \a% ... a r ] . 
A similar parametrization can be found for Toeplitz matrices. 
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The Cramer-Rao bound is then given by 

CRB S (X) =tr(A ( ! t J- 1 (a)A c 



(6) 



where A a = S T and 3(a) = - E [V 2 a p(y; a)] is the 

Fisher information matrix of a ifTTl . 

3(a) = A„A T C- 1 AA ( ! t . 

2) Positive semidefinite matrix: Positive semidefinite matri- 
ces can be parameterized as X = MM T , where M <E R" xr . 
Let a = vec(M), so that X = g(a) then 3(a) has the same 
form as above, with A„ = 2sL2i xhe gradient can be writ- 
ten compactly as = [(I n ® M)T + (M ® I„)] where 

M = mat„. r (o:) and T is the matrix form of the transpose 
operator, Tvec(M) — vec(M T ). Note that parametrization by 
M is unique only up to an orthonormal transformation, hence 
3(a) is not invertible in general. The Cramer-Rao bound is 
then given by 



CRB+(X) =tr A^jt(a)A 



(7) 



provided that 



A; = A;j(a)jt(a) 



holds, or equivalently that P}A a = 0, where Pj is the 
orthogonal projection onto the nullspace of 3(a) |fl~). 

For sake of brevity the CRB derivations for Hankel and 
p.s.d. matrices are given in a supplementary note |[T8l . 

V. Experiments and Results 
A. Experiment setup and performance measure 

In the following we consider real-valued matrices of di- 
mension n = p = 100, with Hankel and p.s.d. struc- 
ture respectively. For Hankel structure we generate X ran- 
domly by first creating a matrix with elements from an i.i.d. 
Af(0, 1) and fitting a using Prony's method [19|. Then set 
X = mat„. p (Sg(a)). For p.s.d. structure, X is generated by 
X = MM T , where the elements of M are generated by 
i.i.d. Af(0, 1). We let parameter A = r / min(n,p) 6 (0,1] 
determine the rank, which controls the degrees of freedom. 

The linear sensing operator, modeled by A in 0, is selected 
randomly by ~ A^(0, ^) ifTJl . The measurement noise is 



generated as n~ A/"(0, er 2 I m ). In the experiments, two signal 
parameters are varied: 

1) Signal to measurement noise ratio, 

E[||X 



SMNR 



El 



m 

in 



(8) 



2) Measurement fraction, p £ (0, 1], so that m = \pnp]. 
The signal-to-reconstruction error ratio, or inverse of the 
normalized mean square error (NMSE), 

Enixn 2 P i i 



SRER = 



E 



X|| 2 F 



NMSE' 



(9) 



is used as the performance measure as it increases with SMNR. 
For an unbiased estimator, SRER < EfUXH 2 ,]/ E[CRB(X)]. 
When SRER = dB there is no reconstruction gain over the 
zero solution X = 0. 



B. Results 

The experiments were repeated for 500 Monte Carlo sim- 
ulations. For each run a new realization of X, y and A 
was generated and the Cramer-Rao bounds were computed 
correspondingly. The algorithm was set to terminate when the 
measurement residual J(L, R) stops decreasing. 

Figure [T] shows the performance of ALS when varying 
SMNR for Hankel and p.s.d. matrices, respectively. The mea- 
surement fraction was fixed to p = 0.3 and the relative 
rank was A — 0.03. The algorithm is tested both with and 
without prior knowledge of matrix structure. Without such 
information it quickly approaches the bound for unstructured 
matrices (O. The reconstruction gain is significantly raised 
when exploiting the matrix structure. ALS remains at a SRER 
gap from the bounds © and (0 because it relaxes the problem 
by alternating projections. The gaps are 2.75 and 0.77 dB 
for Hankel and p.s.d., respectively, at SMNR = 10 dB. The 
estimator inserts a proportionally larger bias into the MSE 
compared to the estimator without prior information. 

In Figure |2] the experiment is repeated for varying p, while 
A = 0.03 and SMNR is fixed to 10 dB. As the measurement 
fraction increases the SRER of the algorithm without prior 
matrix structure approaches the corresponding CRB. At very 
low p, below 0.1, the bound (0) tends to break down as the 
rank constraint is violated. Given prior structural information 
the algorithm performs similar to the SMNR case. 

Reproducible research: The MATLAB code to 
reproduce the figures in this paper are available at 
http://sites.google.com/site/saikatchatt/, along with the 
figures for varying A that have been omitted here due to lack 
of space. 

VI. Conclusions 

We have developed an algorithm based on least-squares 
estimation for reconstruction of low-rank matrices in a general 
underdetermined setup. Furthermore it is capable of exploit- 
ing structures, in particular linearly structured matrices and 
positive semidefinite matrices, leading to better performance. 
Through simulations, we found that the algorithm provides 
good performance compared to the Cramer-Rao bound. 
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